//#include <stdio.h>
#include <iostream>
#include <fstream>
#include <mkl_lapacke.h>

#define index_h(i,j) ((i)*(iter_max)+j)
#define index_kr(i,j) ((i)+(j)*M)
//#define N 3
//#define lda 3
extern "C" void orth_vec(float *P, float *rR, int m, int n,  int tid); 

extern "C" void qr(float *A, int M, int N, int iter_max, float *Q, float * R){
    //N=kR;
    // 输入矩阵 A
   // float A[M][N] = {{1, 2, 3}, {4, 5, 6}, {7, 8, 10},{0,0,1}};
    float *tau=new float[M*N];
    float *A1=new float[M*N];
    //FILE *fp=fopen("output.txt","w");

    // 声明变量
    int i, j, info, lda=N;
    for (i = 0; i < M; ++i) {
                for (j = 0; j < N; ++j) {
                    A1[i*N+j]=A[index_h(i,j)];
                }
            }
    /*
    fprintf(fp,"\nA=\n");
    for (i = 0; i < M; ++i) {
        for (j = 0; j < N; ++j) {
            fprintf(fp,"%5.4e ",A1[i*N+j]);
        }
        fprintf(fp,"\n");
    }*/
    // 执行 QR 分解
    info = LAPACKE_sgeqrf(LAPACK_ROW_MAJOR, M, N, A1, lda, tau);

    if (info == 0) {
//        printf("上三角矩阵 R：\n");
        for (i = 0; i < M; ++i) {
            for (j = 0; j < N; ++j) {
                if (j < i) {
                    R[index_kr(i,j)]=0.0;
                } else {
                    R[index_kr(i,j)]=A1[i*N+j];
                }
            }
        }
//    fprintf(fp,"\nR=\n");
    for (i = 0; i < M; ++i) {
        for (j = 0; j < N; ++j) {
  //          fprintf(fp,"%5.4e ",A1[i*N+j]);
        }
    //    fprintf(fp,"\n");
    }        info = LAPACKE_sorgqr(LAPACK_ROW_MAJOR, M, N, N, A1, N, tau);

      //  printf("正交矩阵 Q：\n");
        if (info == 0) {
            for (i = 0; i < M; ++i) {
                for (j = 0; j < N; ++j) {
                    Q[index_kr(i,j)]=A1[i*N+j];
                }
            }
        } else {
            printf("QR 分解失败！\n");
            return ;
        }
    } else {
        printf("QR 分解失败！\n");
        return ;
    }
    delete[] tau;
    delete[] A1;
    

    Q[index_kr(0,N)]=1.0;
    for(i=1;i<M;i++) Q[index_kr(i,N)]=0.0;
    orth_vec(Q, Q+N*M, M,  N, 0); 
    /*
    fprintf(fp,"\nQ1=\n");
    for (i = 0; i < M; ++i) {
        for (j = 0; j < N; ++j) {
            fprintf(fp,"%5.4e ",A1[i*N+j]);
        }
        fprintf(fp,"\n");
    }
    fprintf(fp,"\nQ2=\n");
    for (i = 0; i < M; ++i) {
        for (j = 0; j <= N; ++j) {
            fprintf(fp,"%5.4e ",Q[i+j*M]);
        }
        fprintf(fp,"\n");
    }
    fclose(fp);
*/
    /*
        printf("Upper矩阵 R：\n");
            for (i = 0; i < M; ++i) {
                for (j = 0; j < N; ++j) {
                    printf("%-8g ", R[index_kr(i,j)]);
                }
                printf("\n");
            }
        printf("正交矩阵 Q：\n");
            for (i = 0; i < M; ++i) {
                for (j = 0; j < N; ++j) {
                    printf("%-8g ", Q[index_kr(i,j)]);
                }
                printf("\n");
            }
            */
    return ;
}

